pic = imread('lena_s.png');%'Mars-1.jpg''circle.png'
sigma = 1.6;
sz = size(pic);
g=zeros(sz(1), sz(2), 3, 5);
hog = zeros(sz(1), sz(2), 3, 4);
imshow(pic);
for i = 1:5
    g(:,:,:,i) = imfilter(pic, fspecial('gaussian',[16,16],sigma*(1+2^((i-1)/2))),'replicate');
    if i > 1
        hog(:,:,:,i-1)=g(:,:,:,i)-g(:,:,:,i-1);
    end
end
out1 = zeros(sz(1), sz(2));
dir = [-1,-1; -1,0; -1,1; 0,-1; 0,1; 1,-1; 1,0; 1,1];
for x = 2 : sz(1)-1
    for y = 2 : sz(2)-1
        flag = 0;
        % min
        for i = 1:8
            if hog(x,y,1,2) < hog(x+dir(i,1), y+dir(i,2),1, 1) && hog(x,y,1,2) < hog(x+dir(i,1), y+dir(i,2),1, 3) && hog(x,y,1,2) < hog(x+dir(i,1), y+dir(i,2),1, 2)
                continue;
            else
                flag = 1;
                break;
            end
        end
        if flag == 1
            flag = 0;
            for i = 1:8
                if hog(x,y,1,2) > hog(x+dir(i,1), y+dir(i,2),1, 1) && hog(x,y,1,2) > hog(x+dir(i,1), y+dir(i,2),1, 3) && hog(x,y,1,2) > hog(x+dir(i,1), y+dir(i,2),1, 2)
                    continue;
                else
                    flag = 1;
                    break;
                end
            end
            if flag == 0 && hog(x,y,1,2) > hog(x,y,1,1) && hog(x,y,1,2) > hog(x,y,1,3)
                out1(x,y) = 255;
            end
        else
            if hog(x,y,1,2) < hog(x,y,1,1) && hog(x,y,1,2) < hog(x,y,1,3)
                out1(x,y) = 255;
            end
        end
    end
end
out2 = zeros(sz(1), sz(2));
for x = 2 : sz(1)-1
    for y = 2 : sz(2)-1
        flag = 0;
        % min
        for i = 1:8
            if hog(x,y,1,3) < hog(x+dir(i,1), y+dir(i,2),1, 2) && hog(x,y,1,3) < hog(x+dir(i,1), y+dir(i,2),1, 4) && hog(x,y,1,3) < hog(x+dir(i,1), y+dir(i,2),1, 3)
                continue;
            else
                flag = 1;
                break;
            end
        end
        if flag == 1
            flag = 0;
            for i = 1:8
                if hog(x,y,1,3) > hog(x+dir(i,1), y+dir(i,2),1, 2) && hog(x,y,1,3) > hog(x+dir(i,1), y+dir(i,2),1, 4) && hog(x,y,1,3) > hog(x+dir(i,1), y+dir(i,2),1, 3)
                    continue;
                else
                    flag = 1;
                    break;
                end
            end
            if flag == 0 && hog(x,y,1,3) > hog(x,y,1,2) && hog(x,y,1,3) > hog(x,y,1,4)
                out2(x,y) = 255;
            end
        else
            if hog(x,y,1,3) < hog(x,y,1,2) && hog(x,y,1,3) < hog(x,y,1,4)
                out2(x,y) = 255;
            end
        end
    end
end

for i = 1:5
    
end
% for i = 1:4
%     figure(i);
%     imshow(hog(:,:,:,i)/.255);
% end
figure(1);imshow(pic);
figure(2);imshow(out2-out1);